Resummed Green-Kubo relations for a fluctuating fluid-particle model 
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A recently introduced stochastic model for fluid flow can be made Galilean invariant by intro- 
ducing a random shift of the computational grid before collisions. This grid shifting procedure 
accelerates momentum transfer between cells and leads to a collisional contribution to transport 
coefhcients. By resumming the Green-Kubo relations derived in a previous paper, it is shown that 
this collisional contribution to the transport coefficients can be determined exactly. The resummed 
Green-Kubo relations also show that there are no mixed kinetic-collisional contributions to the trans- 
port coefficients. The leading correlation corrections to the transport coefficients are discussed, and 
explicit expressions for the transport coefficients are presented and compared with simulation data. 



PACS numbers: 47.11. 05.40.-a, 02.70.Ns 

A recently introduced stochastic model for fluid flow 
P, with efficient multi-particle interactions — which 
we will call Stochastic Rotation Dynamics (SRD) — is a 
promising tool for the coarse-grained modeling of a fluc- 
tuating solvent, particularly for colloidal 3] and polymer 
suspensions 0, H, • SRD can be thought of as a "hy- 
drodynamic heat bath" , the details of which are not fully 
resolved, but which provides the correct hydrodynamic 
interactions between embedded particles. In addition to 
its numerical advantages, its simplicity makes it possible 
to obtain analytic expressions for the transport coeffi- 
cients which are valid for both large and small mean free 
paths, something which is very difficult to do for other 
mesoscale particle-based methods. 

In its original form [l|, Q , the SRD algorithm was not 
Galilean invariant at low temperatures, where the mean 
free path. A, is smaller than the cell size a. However, 
as was shown in Refs. Q and Q, Galilean invariance 
can be restored by introducing a random shift of the 
computational grid before every multi-particle interac- 
tion. A discrete-time projection operator technique was 
then used 8] to derive the Green-Kubo (GK) relations 
for the model's transport coefhcients. Using these results, 
explicit expressions for the transport coefficients were de- 
rived in an accompanying paper Q. In particular, it was 
shown that the grid shifting procedure accelerates mo- 
mentum transfer between cells and leads to a collisional 
contribution to the transport coefficients. However, the 
resulting expressions, while accurate, were only approx- 
imate, since it was not possible to sum-up in any con- 
trolled fashion all the terms in the GK relations. Sub- 
sequently, Kikuchi et al used a non-equilibrium ap- 
proach to derive expressions for the shear viscosity which 
differed slightly from those derived in 0. Furthermore, 
while their approach yielded only two — pure kinetic and 
collision — contributions to the viscosity, the analysis of 
the GK formalism presented in Refs. Q and @ suggested 
that there are additional "mixed" contributions. These 



discrepancies led us to re-examine the GK approach. 

In this paper we show that it is possible to resum the 
time series in the GK relation in such a way as to elimi- 
nate all dependence on the particles' space-fixed cell co- 
ordinates. This leads to dramatic simplifications and al- 
lows the exact evaluation of the collisional contribution to 
the transport coefficients. Furthermore, it is shown that 
there are only pure kinetic and collision contributions to 
the transport coefficients, with no cross terms. 

In the SRD algorithm, the fluid is modeled by parti- 
cles with continuous spatial coordinates (t) and veloci- 
ties Vi (t) . The system is coarse-grained into the cells of a 
regular lattice with no restriction on the number of parti- 
cles in a cell. The evolution of the system consists of two 
steps: streaming and collision. In the streaming step, the 
coordinate of each particle is incremented by its displace- 
ment during the time step, r. Gollisions are modeled by a 
simultaneous stochastic rotation of the relative velocities 
(relative to the mean velocity of the particles in a cell) of 
every particle in each cell. As discussed in Refs. and 
m , a random shift of the particle coordinates before the 
collision step is required to ensure Galilean invariance. 
All particles are shifted by the same random vector with 
components in the interval [—a/2, a/2] before the colli- 
sion step. There is a great deal of freedom in how the 
rotation step is implemented, and any stochastic rotation 
matrix consistent with detailed balance can be used. In 
two dimensions, the stochastic rotation matrix is typi- 
cally taken to be a rotation by an angle ±a, with prob- 
ability 1/2 Q. In three dimensions, two collision rules, 
denoted by Models A and B in Ref. 0, have been con- 
sidered. In Model A one performs rotations by an 
angle a about a randomly chosen axis. In Model B |[TTI |. 
rotations are performed about one of three orthogonal 
rotation axes of a cartesian coordinate system. At each 
collision step, one of these three axes is chosen at ran- 
dom, and a rotation by an angle ±a is then performed, 
where the sign is chosen at random. 
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FIG. 1: G{nT) as a function of time step for A/a = 0.01 (o) 
and A/a = 1.0 (□). The bullets (•) are the result Gcinr) 
given in Eq. ||^. Averages were taken 100,000 iterations and 
5 different random number seeds. Parameters: L/a — 64, 
M = 5, and a = 90°. 




FIG. 2: SG(2t) as a function of A/a for rotation angles a = 
60° (o), 90° (□), and 120° (►). The asterisks (*) are the 
theoretical values values for A/a — 0. Time averages over 
400,000 iterations for 5 different random number seeds were 
used to obtain the data. Parameters: L/a = 64 and M = 5. 



Explicit GK relations for the transport coefficients of 
the SRD algorithm were derived in Ref. In particular, 
it was shown that the shear viscosity, j/, is given by [9llll| 

oo ' 

^ = NksT ^ {<^xy{0)a^y{nT)) , (1) 

n— 

where 

o-xyinr) = 

-lY.^^rAnr)A^,y{nT) + Av,x{nT)Aejy{nT)] , (2) 
j 

with A^^. (nr) = ([n + 1] r) - (nr), (nr) = 
^^■([n+l]r) - ^j([n + l]T), and Av^j (riT) 
Vxj{[n + 1]T) — Vxj{nT). ^j{nT) is the cell coordi- 
nate of particle j at time nr, while is it's cell 
coordinate in the (stochastically) shifted frame. The 
prime on the sum indicates that the t = term has the 
relative weight 1/2. The sum in Eq. Q runs over all 
N particles of the system. Here and in the following we 
have set the particle mass equal to one. 

The straightforward evaluation of the GK relations 
presented in Ref. Q leads to three contributions to the 
transport coefficients, which were called the kinetic, ro- 
tational and mixed terms. For large mean free path, 
A = Ty/ksT 3> a, the assumption of molecular chaos is 
valid, and the kinetic contribution could be determined 
explicitly. For mean free paths smaller than the cell size, 
however, there are finite cell size corrections, and it was 
not possible to sum these contributions in a controlled 
fashion. The origin of the problem was the explicit ap- 
pearance of the cell coordinate in the stress correla- 
tion functions. 

In fact, the appearance of A^ is troubling, since one 
would not expect this to be the case if the cell shifting 



procedure really does restore Galilean invariance. The 
key to resolving this dilemma is to realize that a proper 
resummation of the GK relations removes this depen- 
dence. In particular, by canceling ^-dependent terms in 
successive contributions to the time series in Eq. ^ and 
using stationarity |12| . it can be shown that transport 
coefficients are given by the same GK relations, but with 
the stress tensor axyinr) = a^'^jP{nT) + (f^°*{nT), with 

j 

and 

^lyinr) = XI Bjy{nT)vjx{nT), (4) 
j 

where Bjpinr) = ^^"^ ([n-f 1]t) - C|^(nT) - TVjpinr). 
Note that the new stress tensor does not depend on ^, 
the space-fixed cell coordinates of the particles. It can 
be shown 12], and has been verified numerically, that 
{Bia) = and that all correlations of the _B- fields with 
the particle velocities in the stress correlation fimctions 
factorizc. Furthermore, 

{Bia{nT)Bji}{mT)) = 

^0/3(1 + ^ij) [2(5n,m — 5n,m+l — ^n^m-l] , (5) 

so that the B's are uncorrelated for time lags greater than 
one time step. These relations imply that there are only 
two — a pure kinetic and a pure rotational — contributions 
to the transport coefficients. Relation ||SJ) is of central im- 
portance, because it contains all the geometrical features 
of the grid that contribute to the transport coefficients, 
and is independant of specific collision rules and particle 
properties. 
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the number of particles in a cell one finds 
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FIG. 3: Simulation results for the normalized rotational con- 
tribution to a) the kinematic viscosity, VrotTja^ , and b) the 
thermal diffusivity, DT,rotT /a? , as a function of the collision 
angle a. The bullets (•) are the diagonal, the squares (■) 
the off-diagonal, and the triangles (a) the total contribu- 
tion to the rotational viscosity and thermal diffusivity. The 
solid lines are the theoretical predictions. The data were ob- 
tained by time averaging over 360,000 iterations. Parameters: 
L/a = f6. A/a = O.f , M = 3 and r = f . 



Using these results in (Q, the viscosity can be written 
a.s V ^ Vkin + Vrot-, with 
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Assuming molecular chaos, it is straightforward to eval- 
uate the kinetic contribution to the shear viscosity. If, 
in addition, it is assumed that the number of particles in 
any cell is Poisson distributed at each time step, with an 
average number M of particles per cell, and average over 



2D 
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(8) 



.(M- l + e-A^)sin^(a) 

in two dimensions, where 

Gcinr) ^ {a%-{Q)a%-{nT))clN{kBTf 

= [l~2sin"(a)(Af- 1-f e-^'O/A^]"- (9) 

The index C indicates that molecular chaos was assumed 
when performing the averages. The corresponding result 
for Model A in three dimensions is 
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{M -1 + e-^0[2 - cos(a) - cos(2a)] 
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(see also Refs. [Til El)- 

Eqs. (|HJl and H1U|I are the same results one would 
obtain in the Chapman-Enskog approximation 0. For 
small mean free path, however, there are significant con- 
tributions to Vkin which are neglected in this approx- 
imation. They arise from correlations between parti- 
cles which are in the same (shifted) cell at more than 
one time step. Fig. ^ contains a plot of the G{nT) = 
{a^'-y {<d)d^™ {tit)) / N {ksTY , in two dimensions for a = 
90°. The bullets are the result Gcinr) given in 0, and 
the open squares are simulation data for A/a = 1; the 
agreement shows that for this value of the mean free path, 
Eq. (jSJ provides an excellent approximation for Vkin ■ 

On 

the other hand, data obtained for A/a = 0.01 (o) exhibit 
much larger correlations for nr > 2. Fig. |21 contains a 
plot of the relative difference, SG{2t) = G(2t) - G'c(2t) 
as a function of the mean free path. While it is rather 
difficult to evaluate these corrections analytically for gen- 
eral A, we have calculated 5G{2t) in the A ^ limit. The 
results of this calculation, which are shown in Fig. [3 are 
in excellent agreement with the numerical results. 

There are corrections of this type at small A/a for all 
the transport coefficients, and it is important to note 
that they provide a particularly large contribution to the 
bare self-diffusion coefficient [aiBll3- The effect of these 
correlations on the value of the viscosity are less signif- 
icant and only visible at intermediate mean free path, 
since they vanish for large A and are small compared to 
the dominant coUisional contribution for A <C a. For 
A/a — 0.4, the correlations at ti = 2 make an additional 
contribution of approximately 12% to the total viscosity. 

The rotational contribution to the viscosity is easy to 
evaluate, since, as can be seen from Eq. O, only stress 
correlation functions at equal time and for a time lag 
of one time step are required. Another simplifying fea- 
ture is that because of momentum conservation, the di- 
agonal (from i = j) and off-diagonal (from i ^ j) con- 
tributions to lyrot in O obey the relation i/. 
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Using this result and relation (O, and 



avera ging over the number of particles in a cell, one ob- 
tains [la 



6dT \ 
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M 



[1 — cos(a)] 



(11) 



for all the collision models we considered (the standard 
model in d = 2 and both models A and B in 3d ^3)- ^'i- 
agrees with the result of Kikuchi et al obtained 
using a different non-equilibrium approach in shear flow, 
but deviates slightly for small M from the result given in 
Refs. and jllj. Result Hll|) is compared with simula- 
tion data for the rotational contributation to the viscosity 
in Fig. Ot. 



The GK relation for the thermal diffusivity, Dt, de- 
rived in Ref. 0, 0| can be resummed in a similar fash- 
ion. In particular, it can then be shown that Dt = 
DT,kin + DT.rnt- DT,kin was Calculated in 2d in Ref. 
and in 3d in [ll| neglecting fluctuations in the number 
of particles in a cell. As for the viscosity, it is straight- 
forard to include particle number fluctuations by aver- 
aging the contributions to the heat-flux correlation func- 
tions over the number of particles in a cell; the result- 
ing expression will be given elsewhere [T^ . The relation 

which follows from energy 
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conservation, can be used to show that the rotational 
contribution to the thermal diffusivity is 
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1 - ^ ) [1 - cos(a)] 



(12) 
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to leading order for large M . Note that in contrast to the 
viscosity, the rotational contribution to the thermal dif- 
fusivity is 0(1/M), so that the corrections to Dt at small 
A/a arising from correlated coUisons are more important 
than for the viscosity. Simulation results for Dt^toi are 
compared with Eq. (fT^ in Fig. ^p. 

It is now clear that the random shift procedure in- 
troduced in Refs. and 8] not only restores Galilean 
invariance, but also enables an exact evaluation of the 
collisional contribution to the transport coefficients and 
clarifies several aspects of the underlying algorithm. In 
addition, the current approach justifies in detail several 
assumptions used in the non-equilibrium calculation of 
Kikuchi et al UOl which led them to the same, correct 



results for the shear viscosity derived here using GK re- 
lations. An advantage of the current approach is that 
it can be used to analyze the transport coefficients of 
the longitudinal modes, namely the bulk viscosity and 
thermal diffusivity, which are hard to calculate in a non- 
equilibrium approach [13| . It can also be used to show 
that the bulk viscosity is equal to zero 0, 0| . 
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